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PRECONDITIONED  ITERATIVE  METHODS  FOR  NONSELFAOJOINT 
OR  INDEFINITE  ELLIPTIC  BOUNDARY  VALUE  PROBLEMS 


Janes  H.  Bramble 
Mathematics  Department 
Cornell  University 
Ithaca,  New  York 
U.S.A. 


and 


Joseph  E.  Pasclak 
Brookhaven  National  Laboratory 
Upton,  New  York 
U.S.A. 

a. u.  ^ r  F 

consider  a  Gal erkin- Finite  Element  approximation 
to  a-general  linear  elliptic  boundary  value  oroblem 
which  may  be  nonselfadjoint  or  Indefinite.  wT^e/ 
show  how  to  precondition  the  equations  so  that  the 
resulting  systems  of  linear  algebraic  equations 
lead  to  Iteration  procedures  whose  Iterative 
convergence  rates  are  independent  of  the  number  of 
unknowns  In  the  solution. 


1.  INTRODUCTION. 

In  recent  years,  the  application  of  Iterative  methods  to 
preconditioned  linear  systems  has  been  extremely  successful  In  a  variety 
of  complex  physical  applications  [3,1$].  Many  articles  are  available 
In  the  literature  which  report  on  the  favorable  performance  of  such 
methods  [34*10,12]. 

The  two  aspects  of  a  resulting  algorithm  consist  of  the 
preconditioner  and  the  underlying  Iterative  method  [1,8,12].  Various 
Iterative  methods,  the  most  popular  being  the  conjugate  gradient  (CS)  and 
certain  normal  forma  of  the  C6  method,  have  been  considered  extensively 
both  from  a  theoretical  and  an  experimental  viewpoint  (see  [10]  and  the 
refarences  therein).  It  has  been  demons tri ted  that.  In  general, 

same  theoretical  convergence  retee 


Iterative  algorithms  with  the 


converge,  in  practice,  at  about  the  same  rate^.  The  question  of  choosing 
an  appropriate  preurnditioner  Is  much  more  difficult.  The  - 
preconditioner  must  In  some  way  be  similar  to  the  inverse  of  the  system 
which  Is  being  solved.  Consequently,  the  evaluation  of  the 
preconditioner  usually  requires  the  solution  of  a  system  of  equations  and 
so  if  the  method  Is  to  result  in  an  improvement  of  computational 
efficiency,  the  preconditioner  must  have  some  property  which  makes  It 
easier  to  solve  than  the  original  system.  The  iterative  convergence 
rate  of  the  algorithm  is  extremely  sensitive  to  the  choice  of 
preconditioner.  Indeed,  the  choice  of  a  more  anproorlate 
preconditioner  may  reduce  the  number  of  iterations  by  an  order  of 
magnitude  or  more  in  a  given  problem. 

In  this  paper  we  Illustrate  some  techniques  for  analysing 
preconditioned  Iterative  methods  for  nonsymmetric  problems.  We  will 
discuss  the  problem  of  choosing  an  appropriate  preconditioner  and  study 
two  different  Iterative  algorithms.  Typical  finite  element 
discretization  of  an  elliptic  boundary  value  problem  leads  to  a  matrix 
problem 

(1.1)  Me  »  d. 

where  M  Is  the  "stiffness"  matrix  associated  with  the  discretization 
and  Is  nonsingular  and  c  and  d  are  vectors.  We  seek  a 
precondl tl oner  Mj1  such  that  Mj  is  symmetric  positive  definite, 

(Mj)"1  is  easier  to  compute  than  (M)*\  and  •approximates 

In  some  sense"  (M)’*.  System  (1.1)  can  of  course  be  replaced  by  the 
equivalent  system 

(1.2)  M*  M'1  Me  •  N*  Mj1  If]  d  . 

Tht  matrix  M'  3  MC  MT*  N  1,  symmetric  positive  definite  and  the  . 
first  algorithm  Is  defined  by  applying  the  conjugate  gradient  method 
to  (1.2).  Alternatively,  (1.1)  Is  equivalent  to  the  problem 

(1.3)  Mj1  M*  M^  Me  ■  Mj1  M*  M^1  d. 

^  The  ntmfcer  of  Iterations  to  reach  a  desired  accuracy  may  vary  by  at  most 
a  factor  of  five  [6,10]. 


The  matrix  M"  2  M  although  not  usually  symmetric.  Is  a 

symmetric  operator  with  respect  to  the  Inner  product  defined  by 

«w,v»  ■  (Mj  w)  •  v  . 

The  CG  method  can  be  applied  to  (1.3)  In  the  «*,»»  Inner  product 
and  leads  to  Algorithm  II  of  Section  2.  Our  analysis  suggests  that  the 
preconditioned  Iterative  method  based  on  (1.3)  Is  more  robust  than  that 
based  on  (1.2)  since  results  for  (1.2)  require  additional  hypotheses.  In 
fact,  we  have  not  been  able  to  obtain  results  for  the  scheme  based  on 
(1.2)  unless  the  elements  used  in  the  methods  are  of  "quasi -uni form" 
size. 

We  shall  present  two  general  theorems  which  can  be  used  to  derive 
certain  discrete  stability  estimates.  Such  estimates  lead  to  bounds  on 
the  Iterative  convergence  rates  of  algorithms  for  finding  the  solution  of 
matrix  equations  resulting. from  the  finfte  element  discretization  of 
elliptic  boundary  value  problems  which  may  be  nonsymmetrlc  and/or 
indefinite.  We  show  how  these  general  results  can  be  applied  In  a 
finite  element  approximation  to  the  Poincare  proa i am.  Both  strategies 
depend  upon  a  priori  stability  estimates  for  the  continuous  problem  and 
use  the  approximation  properties  of  the  discretization  to  derive  the 
stability  estimate  for  the  matrix  problems. 

The  first  theorem  leads  to  a  strategy  which  uses  a  positive  definite 
synraetric  problem  as  a  preconditioner  for  a  more  complicated 
nonsymmetrlc  and/or  indefinite  problem.  The  problem  of  the  efficient 
solution  of  positive  definite  problems,  although  not  completely  solved, 
has  been  extensively  researched.  For  example,  matrices  corresponding 
to  positive  definite  symmetric  problems  often  have  certain  diagonal 
dominance  properties  which  Imply  that  various  sparse  matrix  packages 
[9,11]  can  be  used  for  their  solution.  Also,  there  are  "fast  solver" 
algorithms  available  for  certain  elliptic  problems  on  a  variety  of 
domains  [5,14,15].  Our  analytical  results  guarantee  that  the  Iterative 
convergence  rate  for  our  algorithms  Is  Independent  of  the  number  of 
unknowns  In  the  system.  Thus  the  cost  of  convergence  to  a  given 
accuracy  grows  linearly  with  the  size  of  the  problem. 

The  first  strategy  is  applicable  to,  fbr  example,  problems  trftere 
the  differential  operator  A  can  be  decomposed  Into  a  symmetric  positive 
definite  operator  L  and  a  compact  (but  not  small)  perturbation  B.  The 


operators  A,  L,  and  B  are  approximated  by  discrete  operators  A^, 

Lh,  and  derived  by  finite  elements.  The  discrete  approximation 
to  the  solution  u  of  the  original  problem  Is  defined  as  the  solution 
of 

(1.4)  (Lh  +  Bh)U  *  F* 

Problem  (1.4)  can  be  replaced  by  the  equivalent  problem 

<’-5>  m1  <h,  *  »h)u  ■  r  ■ 

We  derive  the  appropriate  stability  estimates  for  (1.5)  which  guarantee 
that  the  CG  method  applied,  with  respect  to  «*,•»,  to  (1.3)  converges 
at  a  rate  Independent  of  the  number  of  unknowns  In  the  discretization.  In 
addition,  the  stability  results  yield  Immediately  estimates  for  the 
discretization  error  u-U. 

We  give  a  second  theorem  which,  under  additional  hypotheses, 
provides  another  stability  estimate.  This  estimate,  under  a  further 
restriction,  can  be  used  to  show  that  the  CG  method  applied  to  (1.2) 
converges  to  the  solution  of  (1.2)  at  a  rate  which  Is  Independent  of  the 
number  of  unknowns  In  the  discretization. 

An  outline  of  the  remainder  of  the  paper  Is  as  follows.  In  Section  4 
we  describe  two  conjugate  gradient  algorithms  for  matrix  problems. 

Section  3  gives  some  preliminaries  and  notation  to  be  used  In  the  paper. 

In  Section  4  we  state  the  type  of  estimates  needed  to  guarantee  rapid 
convergence  for  some  Iterative  methods  for  solving  nonsymmetric  and/or 
Indefinite  problems.  Two  theorems  used  to  derive  the  stability  estimates 
are  given  In  Section  5.  In  Section  6  we  apply  the  theorems  to  a  finite 
element  approximation  of  a  general  elliptic  boundary  value  problem. 

Finally  In  Section  7  we  apply  a  stability  estimate  to  bound  the 
discretization  error. 


2.  CONJUGATE  GRADIENT  ALGORITHMS. 

We  describe  the  algorithms  which  result  from  applying  the  conjugate 
gradient  method  to  the  preconditioned  sy sterna  (1.2)  and  (1.3).  In  either 
cue  we  assume  that  we  are  given  an  Initial  approximation  Cq  to  the 
solution  c  of  (1.1)  and  the  Iterative  algorithm  producu  a  sequence  of 


Iterates  e.|  for  1  >  G.  We  stop  th«  Iterative  procedure  when  the 
residual  error  d-Mc  becomes  sufficiently  snail.  We  note  that  applying 
the  conjugate  gradient  method  to  preconditioned  systens  as 
Illustrated  In  the  following  algorithms  Is  not -novel  however  we  Include 
the  details  for  completeness. 

Applying  the  conjugate  gradient  method  to  (1.2)  gives  the  following 
algorlthai: 

ALGORITHM  I.  M-  •  M*  M^’  M^'  M 

(1)  Opffnt  r0  •  P0  •  M{  Mj1  Mj1  (d-MCj). 


(2)  For  1  >  0  dpflm 


r1  *  h 


°i  *  on 


<=1.1  ■  e1  ♦  “I  >1 


Vi  ■  rt  -  *!  "  "I 
("•  W»Pl 

•l  (M1  p,J*  p, 

pi*i  ■  Vi  -  *1  P1 


Applying  the  conjugate  gradient  method  In  the  «•••»  Inner 
product  to  (1.3)  gives  the  following  algorithm: 

ALGORITHM  II.  r  •  Mj1  M*  MJ1  M  . _ 


(T)  OtflM  Pg  .  Pg  •  M,-’  M*  (if’CMfcj). 
(2)  for  1  >0  define 

<Vl>  •  p, 

■»  * 

*M  *«l  **1»1 
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3.  PRELIMINARIES  AND  NOTATION. 

Throughout  this  paper  we  shall  be  concerned  with  solving  boundary 
value  problems  on  a  bounded  domain  a  contained  In  R2  with 
boundary  r  .  To  state  our  stability  estimates*  we  shall  make  use  of 
various  spaces  of  functions  defined  on  Q  .  The  space  L2(fl)  Is  the 
collection  of  square  Integrable  functions  on  a  ;  that  1$*a  function  f(x) 
defined  for  (x*y)  In  Q  Is  In  L2(Q)  if 

/  f(x*y)2  dxdy  <  •  . 

a 

The  l2(a)  Inner  product  Is  defined  by 

(f.g)  s  /  f(x.y)  g(x*y)dxdy  fbr  f ,  g  e  L2(Q). 

a 

We  shall  also  use  the  Sobolev  space  H^Q).  Loosely,  a  function  f  Is 
In  If  f,  and  are  all  In  L2(Q).  Thus  fbr 

functions  In  we  can  define  the  Dlrlch'let  form  by 

«'•«> » /  JMr  ty**  • 

Me  shall  also  denote  the  l2(r)  Inner  product  by 

<f.9>  i  /  fg  di  . 

r 

Fbr  any  positive  Integer  r,  the  Sobolev  space  of  L2(fl)-ftjnct1ons 
whoso  r**1  order  partial  derivatives  belong  to  l2(a)  will  bo  dsnoted  by 

Mi  also  let  C  and  Cj  ftr  1  >  0  denote  positive  constants.  The 
values  of  C  and  C,  nay  be  different  In  different  places  however  C 
and  C1  shell  always  be  Independent  of  the  mesh  parameter  h  defining 


ths  approximation  method.  Thus  C  and  will  always  be  independent 
of  the  number  of . unknowns  in  the  discretization. 

To  define  the  approximation  of  later  sections  we  shall  need  a 
collection  of  finite  element  approximation  subspaces  {$h>,  0<h£l, 
contained  In  H^(ft).  Typically,  finite  element  approximation  subspaces 
are  defined  by  partitioning  the  domain  a  into  subregions  of  size  h  and 
defining  Sh  to  be  the  set  of  functions  which  are  continuous  on  &  and 
piecewise  polynomial  when  restricted  to  the  subregions  (see  [4,7,17] 
for  details).  For  example,  one  could  partition  ft  into  triangles 
of  size  h  and  define  Sh  to  be  the  functions  which  are  continuous 
on  ft  and  linear  on  each  of  the  triangles.  Alternatively,  ft  could 
be  partitioned  Into  rectangles  and  could  be  defined  to  be  the 
functions  which  are  continuous  on  ft  and  bilinear  on  each  of  the 
rectangles. 

4.  ESTIMATES  FOR  THE  CONJUGATE  GRADIENT  METHOD. 

Our  analysis  of  Iterative  algorithms  for  preconditioned  systems  is 
based  on  stability  estimates  for  the  continuous  or  ncndiscrete  problem 
and  the  error  estimates  between  the  continuous  solutions  and  their 
discrete  approximations.  To  study  the  properties  of  the  solutions  of 
boundary  value  problems  In  partial  differential  equations.  It  is 
natural  to  consider  operators  In  their  basis  free  representations  since 
complete  sets  of  basis  functions  are  usually  too  complex  to  be  of  much 
practical  value.  Consequently,  It  Is  natural  to  think  of  the  process 
of  solving  for  the  discrete  solution  of  the  finite  element  equations  as  a 
basis  free  operator  on  the  finite  element  subspace  $h  of  H'(ft)  .  We 
represent  differential  and  solution  operators  by  the  notation  A,  B, 

L,  or  T  whereas  their  discrete  counterparts  shall  be  respectively 
denoted  A^,  B^,  Lh  and  Th. 

The  C6  method  can  be  applied  to  find  the  solution  *X  of  the  problem 
(4.1)  lh  X  •  Y 

where  Is  a  symmetric  positive  definite  operator  with  respect  to  some 

Inner  product  (cf.  [13]).  The  C6  algorithm  requires  an  Initial  guess  X*  and 

produces  an  approximation  x  to  X  after  n  Iterative  steps.  It  Is' 

n 


well  known  that 

✓  * 

(4.2)  11X-X„1IHI2  . 

where  y  is  the  condition  number  for  and  Is  defined  to  be  thj 
ratio  of  the  largest  eigenvalue  of  Lh  to  the  smallest.  Vie  note  that 
If  Lh  satisfies  the  Inequality 

(4.3)  C0  M*  <  (Ifj  V) ,  W)R  <  C,||M||5  for  >11  M  e  S^. 

where  (•»•)„  denotes  the  H-inner  product*  then  the  condition 
number  y  is  bounded  by  C-j/Cq.  Thus  estimates  of  the  type  (4.3) 
in  conjunction  with  (4.2)  lead  to  convergence  estimates  for  the  CG 
method  applied  to  (4.1). 

The  problem  of  finding  the  finite  element  solution  in  the  examples 
of  later  sections  can  be  reduced  to  solving  for  the  solution  X  of  a 
nonsingular  operator  equation 

(4.4)  \  X  -  Y 

where  ^  Is  a  nonsymmetrlc  and/or  nonpositive  operator  on  $h-  We 
shall  first  precondition  the  system,  multiply  by  the  adjoint  and 
then  apply  the  CG  method  In  the  appropriate  Inner  product. 

We  assume  that  we  have  a  symmetric  positive  definite  operator 
Th  defined  on  Sh  for  a  preconditioner.  The  types  of  preconditioners 
for  which  we  can  get  analytic  results  will  be  described  In  later  sections. 

We  note  that  problem  (4.4)  can  be  replaced  by  the  problem  of 
finding  X  In  S  satisfying 

(4.5)  *h  Th  Th  X  -  Th  Th  Y 

where  A*  Is  the  l2(fl)  -  adjoint  of  The  CG  method  with  respect 
to  the  L2(Q)  Inner  product  can  be  used  to  solve  (4.S).  The 
convergence  rate  of  the  resulting  algorithm  Is  bounded  by  (4.2)  In 
the  l2(S)  norm  where  y  la  bounded  by  Cj/Cg  far  any  Cq  and 
t|  satisfying 


(4.6) 


for  all  W  s  S 


:°“^l2(S!)  ”  l|T"A|'  "llL2(a)lC'll““L2 


(Q) 


In  certain  applications,  estimate  (4.6)  can  be  used  to  derive  bounds  on 
the  iterative  convergence  rate  of  Algorithm  I.* 

Alternatively,  problem  (4.4)  is  also  equivalent  to  the  problem 
of  finding  X  in  satisfying 


(4.7)  Th  Ah  Th  *h  X  “  Th  Ah  Th  Y  ‘ 


The  operator  B  i  A*  A^  is  symmetric  positive  definite  in  the 

inner  product  (T’1  W,  V).  Applying  the  CG  method  to  the  solution  of 

(4.7)  in  this  inner  product  gives  an  algorithm  which  converges  at  a  rate 
described  by  (4.2)  where  y  <  C-j/Cq  for  any  CQ  and  satisfying 

(4.8)  Cgd*1  W,W)  <  (Th  AhW,  A^)  <  C^T*1  W.W)  for  all  W  e  Sh  . 

In  applications,  estimate  (4.8)  is  used  to  derive  iterative  convergence 
rates  for  Algorithm  II. 


5.  STABILITY  THEOREM. 

In  this  section  we  give  general  results  which  can  be  used  to  derive 
estimates  of  the  form  (4.6)  and  (4.8). 

Theorem  1.  Let  R  be  a  continuous  operator  and  Rf1  be  its  discrete 
approximation.  Assume  that  the  following  stability  and  error  estimates 
hold: 


(5.1) 

M  1,  <  e(||(i*Rh)  el,  *  M  2/  . 

h‘(o)  h'(o)  i(n) 

}  for  all  0  e  Sh 

For  any 

c  >  0  there  exists  C£  such  that 

• 

(5.2) 

M  2  <  Cr  2  ♦  1 

V(Q)  -  e  LZ(H)  H '  (Q) 

for  all  #  c  h\q). 

(5.3) 

B(R.RJ#B  9  <  ch  <1  for  all 

h  V(n)  ~  V(n) 

♦  c  HT(n). 

Then  there  exists  hQ  >  0  such  that  for  h  <  hQ 

(5.4)  |jejj  ,  i  C|1(I+R.)0|]. ,  for  all  6  e  S, 

H  (Q)  n  H‘(n)  *  1 

Remark  1 .  Estimate  (5.4)  combined  with 


(5.5)  H(I+Rh)e||  <  diei|  ,  for  all  0eS 
n  H‘(n)  h'(q)  h 

guarantees  a  uniform  (independent  of  h)  iterative  convergence  rate  for 
the  CG  iteration  for  the  solution  of 

(I+Rh)*  (I+Rn)U  »  F 

where  *  denotes  the  adjoint  with  resoect  to  the  (n)  inner 
oroduct.  In  our  finite  element  applications,  I+R^  *  ThAh  and  . 

coM2,  <  (t:1  8.e)  <  C.II6II2,  for  all  6c  S  . 

0  H1  (Q)  "  1  H 1  (IS)  " 

Thus  (5.4)  and  (5.5)  will  imply  (4.8)  for  the  particular  examples  of 
the  next  section. 

12  12 
Theorem  2.  Let  T  and  T  be  continuous  operators  and  Th  and  Th 

be  their  corresponding  discrete  approximations.  Assume  that  the  following 
three  estimates  hold: 


(5.6)  CJT1  y|  2  <  «T2  u||  2  <  C^T1  u|)  ,  for  .11  u  c  l2(n). 

u  i*(n)  L*(n)  1  l‘(si) 

(5.7)  ufl^  i  0,2  1^2  for  a11  u  e  • 

(5.8)  KtJ)-1  U|L2  <  Ch2|U|L2  for  an  U  e  Sh  . 
for  1  •  1,2.  Then 


(5-9>  co«^,n,  *  ,|T'lT^rl  ^(a,-^, 


for  all  U  c  S, 


Remark  2.  Estimate  (5.8)  Is  an  Inverse  property  for  the  operator 
and  in  applications  is  derived  /rom  the  hypothesis  that  the  mesh 
elements  are  of  "quasi  uniform"  size.  Estimate  (5.9)  coincides  with 

(4.6)  when  Ah  •  J 

Remark  3.  The  proofs  of  the  above  two  theorems  are  simple  and 
consequently  will  not  be  included. 


6.  THE  POINCARE  PROBLEM. 

To  illustrate  our  approach  we  consider  a  finite  element 
approximation  of  the  Poincare  problem  in  this  section.  We  consider  the 
following  model  problem: 

3u 


(6.1) 


-Au  +  +  Ku  «  f  in  ft 


!*♦*!¥  +™*°  r 
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where  A  »  — ~  +  — s-  ,  n  and  t  are  respectively  the  normal  and 
3x‘  3y* 

tangential  directions  along  r.  For  simplicity  we  have  considered 
constant  coefficients  in  defining  the  differential  equation  as  well  as 
the  boundary  condition.  Our  results  and  iterative  algorithms  extend  to 
variable  coefficient  problems  without  any  complications.*  We  also  assume 
that  the  solution  of  (6.1)  exists  and  is  unique. 

The  finite  element  approximation  to  (6.1)  can  then  be  defined  by  the 
Galerkln  technique,  ftoltlplylng  (6.1)  by  an  arbitrary  function  *  and 
integrating  by  parts  shows  that  the  solution  u  satisfies 


(6.2) 


D(u,$)  ♦  ♦  K(u,*)  ♦  <6  ♦  ru,*>  -  (/,♦) 


3u 


The  finite  element  approximation  U  to  u  Is  then  defined  to  be  the 
function  U  In  Sh  which  satisfies 


(6.3)  o(u,e)  ♦  (|j,e)  ♦  K(u,e)  ♦  <0  ♦  yU.9>  •  (f.e) 


for  ell  6  c  S 


Equation  (6.3)  can  be  used  to  derive  a  system  of  equations  of  the  form 
(1.1)  defining  the  discrete  solution  U,  i.e.,  using  a  basis *for  $^, 

(6.3)  gives  N  equations  for  the  N  unknowns  defining  U  in  that 
basis. 

To  describe  iterative  methods  for  the  solution  of  (6.3)  and/or 
the  corresponding  matrix  system,  we  shall  need  to  use  some  operator 
notation.  First,  we  consider  the  Neumann  problem 


(6.4) 


w  -  Aw  ■  f  in  8 

3w  n  _ 
3^*0  on  r 


Given  a  function  f  in  l2(8),  the  solution  w  of  (6.4)  is  in  H2(8) 
if  as  we  shall  assume,  r  is  sufficiently  smooth.  We  denote  the 
solution  operator  T  as  the  map  which  takes  f  to  Tf  s  w.  T  is  a 
bounded  map  of  L  (8)  into  H  (8).  The  finite  element  approximation  to 

(6.4)  is  the  function'  W  in  Sh  'satisfying 


(6.5)  0(W,8)  ♦  (W,0)  «  (f,6)  for  all  0  e  Sh  . 


The  discrete  solution  operator  Th  can  then  be  defined  as  the  map  which 
takes  f  to  Th  f  s  W.  Th  Is  a  map  from  l2(a)  onto  Sh  and  the 
following  convergence  estimate  is  well  known  (cf.  [2]): 


(6.6)  Kymi ,  <chifu ,  . 

n  H  (8)  l/(8) 

In  a  similar  manner,  we  can  define  solution  operators  for  the  following 
variational  problems: 

0(X,8)  ♦  (X,4)  ■  (§*■»♦)  +  (k -1)(z.t) 

and 

0(iM)  ♦  (<M)  ■  <8  |jr  ♦  Yw»8>  # 

1  2 

Me  define  the  solution  operators  Rzs  X  end  R  u  s  $.  The  corresponding 


finite  element  approximations  are  given  by  the  solutions  X  and  Y  in 
satisfying 

0{X,8)  +  (X.8)  «  (||,0)  +  (K-l)(z,0)  for' all  8  e  Sh 
and 

D(Y,0)  +  (Y ,0)  *  <6  |£  ♦  vw,0>  for  all  0  e  Sh  . 

respectively.  The  discrete  solution  operators  are  then  defined  by 
1  2 

Rh  z  =  X  and  Rh  u)  =  Y  and  the  following  convergence  estimates  hold: 

(6.7)  krJ.rVj  ,  <  ChW  ,  . 

"  i/(n)  h  (n) 

and 

(6.8)  1(r!-R2HI  2,  ■  <  Ch|lu|  ,  . 

"  l‘(s)  H  (n) 

In  terms  of  operators,  problem  (6.1)  is  equivalent  to 
(I  +  R1  +  r2)u  s  TA  u  »  Tf  . 

The  existence  and  uniqueness  properties  of  solutions  of  <6.1)  can  be  used 
to  show  that  for  any  e  >  0  there  Is  a  constant  C£  such  that 

(6.9)  M  ,  «  C  R(I*fl,+R2)6||  j  ♦  cM  ,  . 

1  l‘(6)  L(Q) 


The  discrete  estimate 


(6.io)  neii ,  <  c{11(i+iJ  +  r?)6|I  ,  **-1011,  >  for  til  •  t  s* 

1  V(a)  -  ^  ^  V(a)  l2(Q)  ^ 

Is  Immediate  from  the  definition  of  .  Problem  (6.3)  can  be  stated 
tn  terms  of  operators  as 

« *  "J  *  #u  *  v* w  •  Ti. f  • 


Applying  Theorem  1  we  get  the  following  stability  estimate: 

✓ 

(6-n)  coKi,q)  i  ita  <-(a,  ,or*" -v 

The  second  Inequality  In  (6.11)  can  be  easily  derived  from  the  definitions. 
The  constants  Cq  and  Cy  In  (6.10)  are  Independent  of  the  mesh 
size  h.  Now  It  Is  easy  to  check  that 

(6.12)  (T*1  W,V)  ■  D(W,V)  ♦  (W,V)  for  all  W,VeSh. 

Comparing  (6.12),  (6.11),  (4.7)  and  (4.8)  implies  that  the  CG  method 
applied  to 

“•13>  VS  TA  “•*  VS  \  r 

converges  with  a  reduction  per  Iteration  which  can  be  bounded 
Independently  of  the  number  of  unknowns. 

Let  PI  and  Nj  respectively  denote  the  "stiffness"  matrices 
corresponding  to  (6.3)  and  (€.6)  In  a  given  basis  A  • 

for  Sy,.  If  the  coefficients  of  a  function  W  In  Sh  In  terms  of  the 
basis  A  are  represented  by  the  vector  c  then 

d  ■  Nj1  A*  Aj1  Ac 

gives  the  coefficients  of  T^  AjJ  T^  U  In  terms  of  A.  Consequently, 
the  sequence  of  vectors  c^  generated  by  Algorithm  II  gives  the 
coefficients  of  the  sequence  of  functions  generated  by  the  C6  method 
appV  d  to  (6.13).  Thus  the  Iterative  convergence  estimates  for 
tf’e  >  ->thod  applied  to  (6.13)  Imply  Iterative  convergence  rates  for 
Al§j  si.  T. 

Th  jve  procedure  Is  an  example  of  mi  Iterative  convergence 
analysis  in  (a) .  No  also  note  that  If  T*  Is  another  discrete 
operator  on  ^  which  Is  spectrally  equivalent  to  T„  In  the  sense  that 

(6.14)  ^(Ty,  «,*)  <  (tj  *J»)  «  C^(Tfc  MtH)  fbr  all  A  «  S* 


•  • 


c 


I 


\ 


then  Tw  can  be  replaced  by  t]  in  (6.11). 

n  n  ?  *  ■  s 

We  next  consider  an  Iterative  analysis  in  L  (n)  based  on 
12  2 

Theorem  2.  Let  T  :  L  (n)  ♦  H  (fl).  denote  the  solution  operator 

for  problem  (6.1)  with  S  «  0,  i.e.,  T^  f  =  u.  The  solution  operator  T1 

satisfies  an  estimate  of  the  form 


(6.15) 


ET1  fi 


<  C,  I1T1  «] 


2  i  M  2  —  U 1  Mi  2 

Lz(n)  “  Lz(a)  1  Lz(n) 


We  have  restricted  to  the  case  of  6  *  0  since  (6.15)  is  well  known  in 
that  case.  Assume  that  both  T^  and  T  can  be  approximated  in  the  same 
finite  element  subspaces  and  let  T^  and  T^  denote  the  corresponding 
discrete  solution  operators.  The  following  convergence  estimates  are  well 
known  for  a  wide  class  of  finite  element  applications  [27]: 


(6.16)  I(t'  -  tVi  ,  <  Cdfyl  ,  . 

n  r(n)  l  vfl) 


We  finally  assume  that  the  Inverse  properties 

(6.17)  E(t!)‘1611  ,  <  Ch'2  ||9||  ,  ,lt  S  . 

h  Lz(n)  Lz(n)  " 

are  also  satisfied.  Estimates  of  the  type  (6.17)  can  usually  be 
derived  from  Inverse  assumptions  for  the  subspaces.  Applying  Theorem  2 
gives  that 

(6.18)  i  HV^'1  ■%)  i  ciMl2(8)  ,or  UcSK  • 

Estimate  (6.18)  guarantees  that  the  C6  method  applied  In  LZ(Q)  for  the 
solution  of 

«••«)  Vhf 

where  A^  ■  (tJ)"*  *  will  converge  to  the  solution  X  at  a  rate  which 
Is  Independent  of  the  number  of  unknowns  In  Sh.  The  resulting  algorithm 
does  not  however  correspond  to  Algorithm  I.  To  guarantee  repld  Iterative 


i 


convergence  rates  for  Algorithm  I  we  must  make  additional  assumptions. 
Again  we  use  the  basis  fi  for  .  If  W  e  we  denote  by  Cy  the 

coefficients  of  W  in  the  basis  6.  We  require  that 

(6.20)  C0(CW#CW)  -  (W*W)L2(n)-  C1(W  fora11  W  e  Sh  * 

Estimate  (6.20)  states  that  the  Gram  or  mass  matrix  is  "equivalent" 
to  the  coordinate  inner  product.  Combining  (6.19)  and  (6.20)  implies 

(6.21)  Cq  |c  I2  <  IM*1  M  cj2  <  C, |c 

for  all  N  dimensional  vectors  c.  Estimate  (6.21)  is  finally  an 
estimate  which  can  be  applied  to  guarantee  uniform  Iterative  convergence 
rates  for  Algorithm  I. 

7.  AN  ESTIMATE  FOR  THE  DISCRETIZATION  ERROR. 

In  order  to  estimate  the  discretization  error  u-U  with  u  and  J 
defined  by  (6.2)  and  (6.3)  respectively,  we  Introduce  the  H1  (£2)- 
projectlon  Ph  onto  S^.  It  Is  defined  for  veH^f,)  by 

(.7.1)  D(Ph  v,9)  +  (Ph  v,0)  »  D(v,0)  +  (v.0),  for  all  9  c  Sh  . 

It  1$  well  known  that  P^  satisfies 

for  v  c  Hr(fl)  and  some  r  >  1  which  depends  on  the  choice  of 
Sh(cf.  [2,7]).  In  view  of  (7.2V  to  estimate  u-U  we  need  only  consider 

P^u-U  .  Hence  we  apply  (5.4)  to  obtain 

|JPhu-U  1|  ,  <  C||(WU(Phu-U)B  ,  , 

"  V(a)“  ^  h  V(a) 

Kith  ll*  nj  ♦  «£  .  Fro*  tin  definition  of  in,  Rj;  „ 

see  that 

(I^HPhM-U)  •  Ph(«1*R2)  (Ph-I)u  • 


Hence 


<  C|jPh(R1+R2)(Ph-I)u|!  , 
n  n  m  '(a) 


from  which  it  follows  immediately  that 
<7-3)  *”*"**'<»  ' 


Thus  using  (7.2)  we  obtain  the  estimate  for  the  discretization  error. 


CU-U|lH’(n, 


<  Ch 


r-1 


M 


Hr(n) 
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